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Abstract 

The partial differential equation of Gaussian diffusion is generalized by us- 
ing the time-fractional derivative of distributed order between and 1, in 
both the Riemann-Liouville (R-L) and the Caputo (C) sense. For a gen- 
eral distribution of time orders we provide the fundamental solution, that 
is still a probability density, in terms of an integral of Laplace type. The 
kernel depends on the type of the assumed fractional derivative except for 
the single order case where the two approaches turn to be equivalent. We 
consider with some detail two cases of order distribution: the double-order 
and the uniformly distributed order. For these cases we exhibit plots of the 
corresponding fundamental solutions and their variance, pointing out the 
remarkable difference between the two approaches for small and large times. 
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1 Introduction 



It is well known that the fundamental solution (or Green function) for the 
Cauchy problem of the linear diffusion equation can be interpreted as a Gaus- 
sian (normal) probability density function (pdf) in space, evolving in time. 
All the moments of this pdf are finite; in particular, its variance is pro- 
portional to the first power of time, a noteworthy property of the standard 
diffusion. 

In this paper we illustrate via fractional calculus two types of generalization 
of this Cauchy problem. One type uses the fractional derivative in the sense of 
Riemann and Liouville (R-L), the other in the sense of Caputo (C). In its uses 
we distinguish between single and distributed orders of fractional derivatives. 
Specifically, we work out how to express their fundamental solutions in terms 
of an integral of Laplace type suitable for a numerical evaluation. Particular 
attention is devoted to the time evolution of the variance for the R-L and C 
cases. It is known that for large times the variance characterizes the type of 
anomalous diffusion. 

The plan of the paper is as follows. 

In Section 2 we write down the two general forms of the time-fractional 
diffusion equation with distributed-order with R-L and C derivatives and the 
Fourier-Laplace representation of the corresponding fundamental solution. 
For this purpose we need to introduce a positive function p(/3) that acts as 
a discrete or continuous distribution of orders. In addition to the particular 
case of a single order /3 with < (3q < 1, we consider two case-studies for 
the fractional diffusion of distributed order: as a discrete distribution we take 
two distinct orders /3 2 with < (3± < (5 2 < 1; as continuous distribution 
we take the uniform density of orders between zero and 1. 

Section 3 is devoted to the time evolution of the variance which is obtained 
from the Fourier-Laplace representation of the corresponding fundamental 
solution, by inverting only the Laplace transform. In the single order case 
we recover the sub-diffusion power-law common to the R-L and C forms; 
for the cases of distributed order we find a remarkable difference between 
the two forms, well visible from their asymptotic expressions for small and 
large times. In section 4 we illustrate our method to get the fundamental 
solutions from their Fourier-Laplace transforms, following the strategy of 
carrying out at first the Fourier inversion and then the Laplace inversion. 
We find instructive to show the graphical representation of the fundamental 
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solutions (in space at fixed times). For the case of fractional diffusion of 
single order, because of the self-similarity of the solutions we limit ourselves 
to show plots of the corresponding solutions at a fixed time t — 1 versus 
x. For the two case-studies of fractional relaxation of distributed order, 
because the self-similarity of the solutions is lost, we provide plots of the 
corresponding solutions versus x at three fixed times, selected as t — 0.1, 
t = 1 and t = 10, contrasting the different evolution for the R-L and the 
C form, in a moderate space-range. We can note how the time evolution of 
the solution in the considered spatial range depends on the different time- 
asymptotic behaviour of the variance for the two forms. 
Finally, concluding remarks are given in Section 5. 

In order to have a self-contained mathematical treatment we have added three 
Appendices: the Appendix A is devoted to the basic notions of fractional 
calculus, whereas Appendices B and C deal with functions of Mittag-Leffler 
and Exponential Integral type, respectively, in view of their relevance for our 
treatment. 



2 The equations for time- fractional diffusion 
of distributed order 

2.1 The R-L and the C forms in space-time domain 

The standard diffusion equation, that in re-scaled non-dimensional variables 
reads 

^u(x,t) = -^u(x,t), xeM, teIR+, (2.1) 

with u(x,t) as the field variable, can be generalized by using the notion of 
fractional derivative of distributed order in tirna^. 

For this purpose we need to consider a function p(f3) that acts as weight for 
the order of differentiation j3 G (0, 1] such that 

p(B) > , and / p(B) d(3 = c> . (2.2) 
J o 

1 We find an earlier idea of fractional derivative of distributed order in time in the 
1969 book by Caputo [5], that was later developed by Caputo himself, see [51 [7] and by 
Bagley & Torvik, see [3] A basic framework for the numerical solution of distributed-order 
differential equations has been recently introduced by Diethelm & Ford [12], Diethelm & 
Luchko [13] andby Hartley & Lorenzo [ [26, 3Tj . 
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The positive constant c can be taken as 1 if we like to assume the normaliza- 
tion condition for the integral. Clearly, some special conditions of regularity 
and behaviour near the boundaries will be required for the weight function 
p(/^jE Such function, that can be referred to as the order density if c = 1, is 
allowed to have 5-components if we are interested in a discrete distribution 
of orders. 

There are two possible forms of generalization depending if we use fractional 
derivatives intended in the R-L or C sense. Correspondingly we obtain the 
time-fractional diffusion equation of distributed order in the two forms: 



9 u(x,t)= f 1 p{(3) t D^ 



dt v w Jo 



d 2 

u(x, t) 



dx 2 



dp, x e E, t > 0, (2.3a) 



and 



r 1 id 2 
/ p(P) t D%uJx,t) dp = —uJx,t) , xGlR,t>0. (2.36) 
JO L J ox z 

From now on we shall restrict our attention on the fundamental solutions of 
Eqs. (2.3a)-(2.3b) so we understand that these equations are subjected to 
the initial condition u(x, + ) = u*(x, + ) = 5(x). Since for distributed order 
the solution depends on the selected form (as we shall show hereafter), we 
now distinguish the two fractional equations and their fundamental solutions 
by decorating in the Caputo case the variable u(x,t) with subscript * as it 
is customary for the notation of the corresponding derivative. 

Diffusion equations of distributed order of both types have been recently 
discussed by several authors: in particular we find the C form e.g. in [7J [9j 
ITUl [TlT l4"6l [51] whereas the R-L form e.g. in [301 EU E2] ■ In some papers the 
authors have referred to the C and R-L forms as to normal and modified forms 
of the time-fractional diffusion equation of distributed order, respectively. 

For a thorough general study of fractional pseudo-differential equations of 
distributed order let us cite the paper by Umarov and Gorenflo [53]. For a 
relationship with the Continuous Random Walk models we may refer to the 
paper by Gorenflo and Mainardi |24j . 



2 For the weight function p(f3) we conveniently require that its primitive P{(3) — 
JqP(P') d/3' vanishes at f3 = and is there continuous from the right, attains the value c 
at [3 = 1 and has at most finitely many (upwards) jump points in the half-open interval 
< (3 < 1, these jump points allowing delta contributions to p((3) (particularly relevant 
for discrete distributions of orders). 
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2.2 The RL and C forms in Fourier-Laplace domain 

The fundamental solutions for the time-fractional diffusion equations (2.3a)- 
(2.3b) can be obtained by applying in sequence the Fourier and Laplace 
transforms to them. We write, for generic functions v(x) and w(t), these 
transforms as follows: 

r+oo 



T {v(x); k} = v(k) 



e^ x v(x) dx, neJR, 



£{w(t);s} 



wis) 



+oo 



-st 



w(t) dt, sGC 



Then, in the Fourier-Laplace domain our Cauchy problems [with u(x, + ) = 
u*(x, + ) = S(x)], after applying formulas for the Laplace transform appro- 
priate to the R-L and C fractional derivatives, see (A. 8') and (A. 9), and 
observing 8(k) = 1, see e.g. [18], appear in the two forms 



SU(K, si 



p(P)s^dp 



1 



) 

p(P)s 1 ~ f3 dp 



p{l3)s p ~ l dp 



u(k, s) 



-K U{K,S 



Then, introducing the relevant functions 



A(s)= / P (P)s l -?dp, 



and 



B(s) = / p(p)s^dp. 



(2.4a) 
(2.46) 

(2.5a) 
(2.56) 



we then get for the R-L and C cases the Fourier-Laplace representation of 
the corresponding fundamental solutions: 

*, . 1 l/A{s) 



U[K,S) 



and 



s + k 2 A(s) k 2 + s/A(s) 
B(s)/s 



(2.6a) 



(2.66) 



k 2 + B(s) 

From Eqs. (2.6a)-(2.6b) we recognize that the passage between the R-L and 
the C form can be carried out by the transformation 



{C: B(s)} 



R-L : 



A(s) 



(2.7) 



We note that in the particular case of time fractional diffusion of single order 
Po (0 < fa < 1) we have p((3) = 5((3 - (5 Q ) hence in (2.5a): A{s) = s 1 "*, 
in (2.5b): B(s) = s*, so that B(s) = s/A(s). Then, Eqs. (2.6a) and (2.6b) 
provide the same result 

u(k, s) = u*(k, s) = — g. (2.8) 

This is consistent with the well-known result according to which the two 
forms are equivalent for the single order case. However, for a generic order 
distribution, the Fourier-Laplace representations (2.6a) (2.6b) are different so 
they produce in the space-time domain different fundamental solutions, that 
however are interrelated in some way in view of the transformation (2.7). 



3 The variance of the fundamental solutions 

3.1 General considerations 

Before trying to get the fundamental solutions in the space-time domain to 
be obtained by a double inversion of the Fourier-Laplace transforms, it is 
worth to outline the expressions of their second moment (the variance) since 
these can be derived from Eqs. (2.6a)-(2.6b) by a single Laplace inversion, 
as it is shown hereafter. We recall that the time evolution of the variance is 
relevant for classifying the type of diffusion. 

Denoting for the two forms 

/+oo f+OO 
x 2 u(x,t)dx, C : al(t) := / x 2 u*(x,t)dx, (3.1) 
-oo J — oo 

we easily recognize that 

R-L : a 2 {t) = -— fi(« = 0,t), C : a 2 (t) = -— = 0, t) . (3.2) 

As a consequence we need to invert only Laplace transforms taking into 
account the behaviour of the Fourier transform for k near zero. 
For the R-L case we get from Eq. (2.6a), 

5(«, a ) = l [l-^ 2 ^ + ...) , 
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so we obtain 

* 2 (s) = -^«(« = 0, s) = . (3.3a) 

For the C-case we get from Eq. (2.6b) 

so we obtain 

^(s) = -^(k = 0,s) = 7 ^. (3.36) 
Except for the single order diffusion, were we recover the well-know result 

a 2 (t) = o*(t) = 2 - , 0</3 <1, (3.4) 

for a generic order distribution, we expect that the time evolution of the 
variance substantially depends on the chosen (R-L or C) form. 

We shall now concentrate our interest to some typical choices for the weight 
function p(/3) that characterizes the time-fractional diffusion equations of 
distributed order (2.3a) and (2.3b). This will allow us to compare the results 
for the R-L form and for the C form. 



3.2 Fractional diffusion of double-order 

First, we consider the choice 

p((3) = Pl 5((3 - A) + p 2 5((3 - fa) , 0< A < & < 1 , (3.5) 

where the constants p\ and p 2 are both positive, conveniently restricted to 
the normalization condition p\ + p 2 = 1. 

Then for the R-L case we have 

A(s) = pi s 1 '^ + p 2 s 1 "* , (3.6a) 
so, in virtue of (3.3a), we have 

a*(s) = 2 Pl s -( 1 +^) + 2p 2 s - {l+M . (3.7a) 
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Finally the Laplace inversion yields, see and compare 



(3.8a) 

Similarly, for the C case we have 

B(s) =p x s Pl +p 2 s 02 , (3.66) 

so, in virtue of (3.3b), 

Finally the Laplace inversion yields, see and compare [PJ, 

2 t& 



Pi \ P2 



P2 T(l + (3 2 ) 
2 



£ — >■ +CX) . 



(3.86) 

Then we see that for the R-L case we have an explicit combination of two 
power laws: the smallest exponent (fix) dominates at small times whereas 
the largest exponent (/? 2 ) dominates at large times. For the C case we have 
a Mittag-Leffler function in two parameters so we have a combination of 
two power laws only asymptotically for small and large times; precisely we 
get a behaviour opposite to the previous one, so the largest exponent (^2) 
dominates at small times whereas the smallest exponent (/3\) dominates at 
large times. 

We can derive the above asymptotic behaviours directly from the Laplace 
transforms (3.7a)-(3.7b) by applying the Tauberian theory for Laplace trans- 
forms^] . In fact for the R-L case we note that for A(s) in (3.6a) s 1- ^ 1 is 
negligibly small in comparison with s 1_ ^ 2 for s — > + and, viceversa, s l ~^ 2 is 
negligibly small in comparison to s 1 ~ /3l for s — > +00. Similarly for the C case 
we note that for B(s) in (3.6b) s^ 2 is negligibly small in comparison to s^for 
s — ► + and, viceversa, s^ 1 is negligibly small in comparison s^for s — > +00. 



3 According to this theory the asymptotic behaviour of a function f(t) near t = 00 and 
t = is (formally) obtained from the asymptotic behaviour of its Laplace transform /(s) 
for s — > + and for s — > +00, respectively. 
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3.3 Fractional diffusion of uniformly distributed order 

Second, we consider the choice 

p(/3) = l, 0<P<1. (3.9) 

For the R-L case we have 



A(s) = s / s' p d/3 



s - 1 
log s 



(3.10a) 



hence, in virtue of (3.3a) 



a 2 (s) = 2 



slog s s 2 \og s 

Then, by inversion, see Appendix C, Eqs (C.15)-(C.16), we get 



(3.11a) 



a 2 (t)=2Kt,0)-z/(M)]~{ 



/2/log(l/t), t^O, 
2t/\ogt , t — > oo 



(3.12a) 



where 



u(t, a) :- 



t a+r 



dr , a > — 1 



/o T(a + r + l) 

denotes a special function just introduced in Appendix C along with its 
Laplace transform. 

For the C case we have 



B(s) 



o 



log s 



(3.106) 



hence, in virtue of (3.3b), 



alts) 



2 logs 



(3.116) 



s s — 1 

Then, by inversion, see Appendix C, Eqs. (C.ll), (C.14), and compare with 
Eqs. (23)-(27) in [9] 



aiit) =2 logt + 7 + e t £ 1 (*) 



2tlog(l/t), t->0, 
21og(t), t^oo 



(3.126) 
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Figure 1: Variance versus t for the uniform order distribution in R-L and C 
forms compared with some single order case. Top: < t < 10 (linear scales); 
Bottom: 10 1 < t < 10 7 (logarithmic scales). 



where 

£i(t) := / du = e / du 

Jt u Jo u + t 
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denotes the exponential integral function recalled in Appendix C, and 7 = 
0.57721... is the so-called Euler-Mascheroni constant. 

For the uniform distribution we find it instructive to compare the time evo- 
lution of the variance for the R-L and C forms with that corresponding to a 
few of single orders. In Fig. 1-Top we consider moderate times (0 < t < 10) 
using linear scales, whereas in Fig. 1-Bottom large times (10 1 < t < 10 7 ) 
using logarithmic scales. 

4 Evaluation of the fundamental solutions 

4.1 The two strategies 

In order to determine the fundamental solutions u(x,t) and u*(x,t) in the 
space-time domain we can follow two alternative strategies related to the 
order in carrying out the Fourier-Laplace in (2.6a) and (2.6b) 

(51) : invert the Fourier transforms getting u(x, s), u*(x, s) and then invert 
the remaining Laplace transforms; 

(52) : invert the Laplace transform getting u(n,t), w*(k, t) and then invert 
the remaining Fourier transform. 

Before considering the general case of time-fractional diffusion of distributed 
order, we prefer to briefly recall the determination of the fundamental so- 
lution u(x,t) (common for both the R-L and C forms) for the single order 
case. 

4.2 The single order diffusion 

For the time-fractional diffusion equation of single order j3 the strategy (SI) 
yields the Laplace transform 

~, s S^ 2 - 1 _| T | S /V2 n , x 

u(x,s) = ^—e W s , 0</3o<l. (4.1) 

Such strategy was adopted by Mainardi [33], [3U [35] to obtain the Green 
function in the form 

u(x,t) =r Po/2 U (\x\/t Po/2 ) , -oo<s<+oo, t>0, (4.2) 

where the variable X := xjt^l 2 acts as similarity variable and the function 
U(x) := u(x, 1) denotes the reduced Green function. Restricting from now 
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on our attention to x > 0, the solution turns out as 

1 1 00 (-x) k 

U (X ) = 5 M ft/2 (x) = -g Hr[ _ w/2+(1 _ ft/2)] 

= — £ ^ r[(ft(* + l)/2] sin[(vr/5 (fc + l)/2] , 
Z7r k=o K - 

where Mp / 2 (x) is an an entire transcendental function (of order 1/(1 — A)/2)) 
of the Wright type, see also [2T1 1221 139] and [37]. 



0.5 



0.4 



0.3 



0.2- 



0.1 - 



P =° 

Po= 1/4 /j 



p = 1/2 



P =3/4 



Pn= 1 





x 



Figure 2: The reduced Green function U(x) = |M / g / 2 (x)) versus x (in the 
interval |x| < 5), for fa = 0, 1/4, 1/2, 3/4, 1. 

Since the fundamental solution has the peculiar property to be self-similar 
it is sufficient to consider the reduced Green function U(x). In Fig. 2 we 
show the graphical representations of U(x) for different orders ranging from 
flo = 0, for which we recover the Laplace density 



U(x) = -e-W, (4.4) 
to fa = 1, for which we recover the Gaussian density (of variance a 2 = 2) 

U(x) = -i= e" x2 /4 . (4.5) 



2,/tt 
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The Strategy (S2): yields the Fourier transform. 

t) = E Po (-kH^) , < p Q < 1 , (4.6) 

where Ep denotes the Mittag-Leffler function, see Appendix B. The strategy 
(S2) has been followed by Gorenflo et al. [20] and by Mainardi et al. [37] 1^2] 
to obtain the Green functions of the more general space-time-fractional diffu- 
sion equations (of single order), and requires to invert the Fourier transform 
by using the machinery of the Mellin convolution and the Mellin-Barnes inte- 
grals. Restricting ourselves here to recall the final results, the reduced Green 
function for the time-fractional diffusion equation now appears, for x > 0, in 
the form: 

\ roo . . \ \ pa+ioo Y(l s) 

U(x) = - J q cos (kx) E Po (-. 2 ) dn = - — l_ ioo T{1 _ M2) * 8 ds , 

(4.7) 

with < a < 1. By solving the Mellin-Barnes integrals using the residue 
theorem, we arrive at the same power series (4.3) of the M- Wright function. 
Both strategies allow us to prove that the Green function is non-negative and 
normalized, so it can be interpreted as a spatial probability density evolving 
in time with the similarity law (4.2). 

For readers' convenience we like to mention other papers dealing with the 
fundamental solutions of fractional diffusion equations (of single ordrer); a 
non exhaustive list includes [D[T51l25ll2Zll2Slll3llSlllSll5nil551and references 
therein. 



4.3 The distributed order diffusion 

Similarly with the single order diffusion, also for the cases of distributed 
order we can follow either strategy (SI) or strategy (S2). In contrast with 
previous papers of our group where we have followed the strategy (S2), see 
[2BI HOI HI] , here we follow the strategy (SI). This choice implies to recall the 
Fourier transform pair (a straightforward exercise in complex analysis based 
on residue theorem and Jordan's lemma) 

C * C e~N rfl/2 , d>0. (4.8) 



d + n 2 2d 1 ' 2 
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In fact we recognize by comparing (4.8) with (2.6a)-(2.6b) that for the RL 
and C forms we have 



. (c = c(s) :=1/A(s) 

\d = d(s):=s/A(s) U " 



c = c(s) := B(s)/s 
d = d(s) := B(s) 



(4.9) 



Now we have to invert the Laplace transforms obtained inserting (4.9) in the 
R.H.S of (4.8). 

For the R-L case we have: 



u(x, s) 



1 



For the C case we have: 



2[sA{s)} 1 / 2 



)S ) = ^M^exp{-N[5( S )]^}. 



cxp 



{-\x\[s/A(s)}^} . 



(4.10a) 



(4.106) 



Following a standard procedure in complex analysis, the Laplace inversion 
requires the integration along the borders of the negative real semi-axis in the 
s-complex cut plain; in fact this semi-axis, defined by s = re ±t7T with r > 
turns out the branch-cut common for the functions s 1 "' 3 (present in A(s)for 
the RL form) and s 13 (present in B(s) for the C form). Then, in virtue of 
the Titchmarsh theorem on Laplace inversion, we get the representations in 
terms of real integrals of Laplace type. 

For the R-L case we get 

u(x,t) = -- J e~ rt Im ju (x,re i7r ) j dr , (4.11a) 

where, in virtue of (4.10a), we must know A(s) along the ray s = re tn with 
r > 0. We write 



A (r c m ^j = p cos(7T7) + ip sin(7T7) , 



where 



p = p (r) = A (r e i7r ) 



7 = l{f) = - arg 

71 



A (re'j 



Similarly for the C case we obtain 

u*(x, t) = — — J e~ rt Im js* (x, re t7T ^j j dr , 



(4.12a) 



(4.13a) 



(4.116) 
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where, in virtue of (4.10b), we must know B(s) along the ray s = re™ with 
r > 0. We write 



B (re CT j = cos(7T7*) + ip* sin(7T7*) 



where 



p* = p*(r) = \B(re™) 
7* = 7*( r ) = - ar § 



7T 



5 re ! 



(4.126) 



(4.136) 



As a consequence we formally write the required fundamental solutions as 



u(x, t) 
and 



J e~ rt P(x,r) dr , P(x, r) = — Im ju (sc, re i7r ) } , (4.14a) 



e r PJx.r) dr , PJx.r) 



-Im{«* (x,re i7r )} , (4.146) 



where the functions P(x,r) and P*(x,r) must be derived by using Eqs. 
(4. 10a)- (4. 14a) and Eqs. (4.10b)-(4.14b), respectively. We recognize that, 
in view of the transformation (2.7), the expressions of P and P* are related 
to each other by the transformation 

P*( r ) r/p(r) , 7„,(r) 1 -7(7-) . ( 4 -15) 

We then limit ourselves to provide the explicit expression for the C form 



P*(x,r) 



1 



lixr 



Im( P y 2 e^7*/2 



exp 



^*l 2 p\l 2 x 



l —p 1 J 2 e -p l J 2xco <' K l*l 2 ) 



2ixr 



sm 



7T7*/2 - pl /2 x sin(7r7*/2) 



(4.16) 

For the R-L form the corresponding expression of P(x; r) is obtained from 
(4.16) by applying the transformation (4.15). 

We note that the fundamental solutions found in this subsection are equiva- 
lent to those obtained by the Authors in [38J by following the strategy (S2) 
after a lengthy manipulation of Mellin-Barnes integrals. 

Herafter we exhibit some plots of the fundamental solutions for the two case 
studies considered in subsection 3.2 in order to point out the remarkable 
difference between the R-L and the C forms. 
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4.4 Plots of the fundamental solutions 



For the case of two orders, we chose {f3\ = 1/4, /9 2 = 1} in order to contrast 
the evolution of the fundamental solution for the R-L and the C forms. 



0.5 



1 1 

R-L 


i i i i i 


i i 

P 1= 1/4 






P 2 =1 - 




t=0.1 






/ t=1 \ 






yy t=io \\ 






i i i i i 







x 



0.5- 




Figure 3: The fundamental solution versus x (in the interval \x\ < 5), for 
the double-order distribution = 1/4, /3 2 = 1} at times t = 0.1, 1, 10. Top: 
R-L form; Bottom: C form. 
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Figure 4: The fundamental solutions versus x (in the interval |x| < 5), for the 
uniform order distribution in R-L and C forms compared with the solutions 
for some cases of single order. Top: t = 1; Bottom: t = 10. 



In Fig. 3 we exhibit the plots of the corresponding solution versus x (in 
the interval |x| < 5), at different times, selected as t = 0.1, t — 1 and 
t = 10. In this limited spatial range we can note how the time evolution of 



17 



the pdf depends on the different time-asymptotic behaviour of the variance, 
for the two forms, as stated in Eqs. (3.12a)-(3.12b), respectively. 

For the uniform distribution, we find it instructive to compare in Fig. 4 
the solutions corresponding to R-L and C forms with the solutions of the 
fractional diffusion of a single order fa = 1/4,3/4, 1 at fixed times, selected 
as t = 1, 10. We have skipped fa = 1/2 for a better view of the plots. 

5 Conclusions 

We have investigated the time fractional diffusion equation with (discretely 
or continuously) distributed order between and 1 in the Riemann-Liouville 
and in the Caputo forms, providing the Fourier-Laplace representation of the 
corresponding fundamental solutions. Except for the case of a single order, 
for which the two forms are equivalent with a self-similar fundamental solu- 
tion, for a general order distribution the equivalence and the self-similarity 
are lost. In particular the asymptotic behaviour of the fundamental solution 
and its variance at small and large times strongly depends on the selected 
approach. We have considered two simple but noteworthy case-studies of 
distributed order, namely the case of a superposition of two different orders 
fa and fa an d the case of a uniform order distribution. In the first case one of 
the orders dominates the time-asymptotics near zero, the other near infinity, 
but f3i and fa change their roles when switching from the R-L form to the C 
form of the time-fractional diffusion. The asymptotics for uniform order den- 
sity is remarkably different, the extreme orders now being (roughly speaking) 
and 1. We now meet super-slow and slightly super-fast time behaviours of 
the variance near zero and near infinity, again with the interchange of be- 
haviours between the R-L and C form. We have clearly pointed out the above 
effects with the figures in sub-section 3.3, in particular the extremely slow 
growth of the variance as t — > oo for the C form. After the analysis of the 
variance, that in practice requires only the inversion of a Laplace transform, 
we have considered the task of the double inversion of the Laplace-Fourier 
representation. For a general order distribution we were able to express the 
fundamental solution in terms of a Laplace integral in time with a kernel 
which depends on space and order distribution in a simple form, see Eqs. 
(4.14)-(4.16). For the two case studies the plots of the fundamental solutions 
(reported in sub-section 4.4) have shown their dependence on the different 
asymptotic behavior of the corresponding variance. 
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Appendix A: The two fractional derivatives 

The purpose of this Appendix is to clarify for the interested reader the main 
differences between the Riemann-Liouville (R-L) fractional derivative and 
the Caputo (C) fractional derivative^ for well-behaved functions f(t) with 
t > 0, exhibiting a finite limit /(0 + ) as t — > + . Denoting with fj, e (0, 1] the 
order, the R-L derivative is defined as 



f - 1 tf-lWi 

jjii f (A . = J r(i - n) dt Jo (t - tY ■ (A . 

' ^ ' dfit) 



dt 

and the C derivative as 



fir) , 

— r / — dr , < ii < 1 , 

t DU{t):={ r ( 1 -^) y ° C*"^ (A2) 



dm 

dt 



li = l. 



The two fractional derivatives are related to the Riemann-Liouville fractional 
integral as follows. The Riemann-Liouville fractional integral is 

tJ " f{t) ]= fJJT) jf/CrMt-r)"- 1 ^, M>0, (A3) 



4 This form is referred to Caputo who used it in the late sixties of the past century, see 
[H [5] . Soon later this derivative was adopted by Caputo and Mainardi in the framework 
of the theory of Linear Viscoelasticity, see [5]. It was mainly with the 1997 CISM chapter 
by Gorenflo and Mainardi [23] and with the 1999 book by Podlubny [47] that such form 
was popularized. 
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(with the convention t J° f(t) = f(t)) and is known to satisfy the semigroup 
property t J^ tJ u — tJ^ +v > with //, v > . For any fi > the Riemann- 
Liouville fractional derivative is defined as the left inverse of the corre- 
sponding fractional integral (like the derivative of any integer order), namely 
t D» t J» fit) = f(t) . Then for /j e (0, 1] we have 

t D» f{t) := t D x t J^ fit) , t D?/(t) := t J 1 ^ t D x fit) , (A4) 

t r t D» fit) = t r tJ 1 ^ tD 1 fit) = tJ 1 t D l fit) = fit) - /(0+) . (A5) 
Recalling the rule 

tDni = t7 -M ^q, 7 >-l, (A6) 

r(7 + 1 - y) 

it turns out for < ji < 1 , 

tDt fit) = t rr [fit) - /(0 + )] = t LT fit) - /(0+) . (A.7) 

Note that for \i — 1 the two types of fractional derivative coincide, the 
constant /(0 + ) playing no role. 

The Caputo fractional derivative represents a sort of regularization in the 
time origin for the Riemann-Liouville fractional derivative and satisfies the 
relevant property of being zero when applied to a constant. 

Let us now consider the behaviour of the above derivatives of non integer 
order with respect to the Laplace transformation^. 

For the Riemann-Liouville derivative of non-integer order /x we have 
£{ t D»fit);s} = s»fis)-giO + ), g(0 + ) :=hm t J l ^f(t), < // < 1. (A8) 
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The Laplace transform of a well-behaved function f(t) is defined as 



f(s)=£{fit);s}:= I e~ s * /(*)<&, *eC. 
Jo 

We recall that under suitable conditions the Laplace transform of the first derivative of 
fit) is given by 

C {tD 1 /(<);«} = a J{a) - /(0+) , /(0+) := lim f(t) . 
1 1 t->o+ 
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For the Caputo derivative of non-integer order fi we have 
C{ t D>tf{ty,s} = s»f(8)- S *- 1 f(0 + ), /(0+):=hm /(*) , < // < 1 . (A9) 

Thus the rule (A. 8) is more cumbersome to be used than (A. 9) since it 
requires initial values concerning an extra function g(t) related to the given 
f(t) through a fractional integral. However, when the limiting value /(0 + ) is 
finite we can see that g(0 + ) is vanishing so the formula (A. 8) simplifies into 

C{ t D"f(t);s} = s»f(s), 0< M <1. (A8') 



For further reading on the theory and applications of fractional calculus we 
recommend to consult in addition to the well-known books by Samko, Kilbas 
& Marichev [19], by Miller & Ross [15], by Podlubny [17], those appeared in 
the last few years, by Kilbas, Srivastava & Trujillo |28j, by Magin [32J, by 
West, Bologna & Grigolini [51], and by Zaslavsky [56J. 



Appendix B: The Mittag-Leffler functions 

B.l The classical Mittag-Leffler function 

Let us recall that the Mittag-Leffler function E^{z) (with fi > 0) is an entire 
transcendental function of order defined in the complex plane by the 
power series 

It was introduced and studied by the Swedish mathematician Mittag-Leffler 
at the beginning of the XX century to provide a noteworthy example of entire 
function that generalizes the exponential (to which it reduces for jj, — 1). For 
details on this function we refer e.g. to [T31 HSl [231 12H1 [361 SZ1 SS] . In particular 
we note that the function E^—x) (x > 0) turns a completely monotonic 
function of x if < \i < 1. This property is still valid if we consider the 
variable x = At M where A is a positive constant. Thus the function E^—Xt^) 
preserves the complete monotonicity of the exponential exp(— At): indeed, 
for < < 1 it is represented in terms of a real Laplace transform (of a 
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real parameter r) of a non-negative function (that we refer to as the spectral 
function) 

E (-Xtn = I /°°e ~ rt Ar "" 1 Sin(/i?r) dr (B 2) 

M ' n Jo A 2 + 2Ar^ cos(/ivr) + r 2 ^ ' v 

We note that as /i — > 1 the spectral function tends to the generalized Dirac 
function <5(r — A). 

We note that the Mittag-Leffler function (B.2) starts at t = as a stretched 
exponential and decreases for t — > oo like a power with exponent — \i: 



l-A— -~exp^- — r , t^0+, 



r(i + M ) r(i + / x) 

r* 4 



I Ar(i-/i) 



(^•3) 

t — > oo . 



The noteworthy results (B.2) and (B.3) can also be derived from the Laplace 
transform pair 

£{E„(-\ns} = -^. (BA) 

In fact it it sufficient to apply the Titchmarsh theorem (s = re 47r ) for deriving 
(B.2) and the Tauberian theory (s — > oo and s — > 0) for deriving (B.3). 

If = 1/2 we have for t > 0: 



E 1/2 (-\Vt) =e Xt erfc(A v / t) ~ l/(\V^t) , t -> oo , (B.5) 
where erfc denotes the complementary error function, see e.g. [2]. 

B.2 The generalized Mittag-Leffler function 

The Mittag-Leffler function in two parameters E^ v (z) (3?{/u} > 0, v 6 C) is 
defined by the power series 

oo k 

It generalizes the classical Mittag-Leffler function to which it reduces for 
v = 1. It is an entire transcendental function of order on which the 
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reader can inform himself by again consulting the references before outlined 
for the classical Mittag-Lefher function 

The function E^ v {—x) (x > 0) is completely monotonic if < fi < 1 and 
v > jj,. Again this property is still valid if we consider the variable x = A 
where A is a positive constant. In this case the asymptotic representations 
as t — > + and t — > +oo read 



1 f 

- A —A r, t^0+, 



(B.7) 



t — ► oo. 



Y(u) T(u + n) 

We point out the Laplace transform pair, see [47] , 

C{t»- 1 E^(-\n,s} = ^^, (B.8) 

with /i, v G 1R + . For < // = z/ < 1 this Laplace transform pair can be used 
to derive the noteworthy identity 

r( 1 -^ w (-Af) = -^£ /1 (-Af) , 0</i<l. (5.9) 

Appendix C: The Exponential integral and 
related functions 

The exponential integral function, that we denote by S\(z), is defined as 

eAz)^ dt= dt. ((7.1) 

Jz t Jl t 

We have used the letter £ instead of E (commonly adopted in the literature) 
in order to avoid confusion with the Mittag-Leffler functions that play a more 
relevant role in fractional calculus. This function exhibits a branch cut along 
the negative real semi-axis and admits the representation 

oo ^ 

= -7-logz- J2 r,|arg2|<vr, (C.2) 

n=1 nn\ 
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where 7 = 0.57721... is the so-called Euler-Mascheroni constant. The power 
series in the R.H.S. is absolutely convergent in all of (D and represents the 
entire function called the modified exponential integral 

Em(z):=/ ——d( = -Y,—, (C.3) 

Thus, in view of (C.2) and (C.3), we write 

S\(z) = — 7 — log 2; + Ein (z) , |argz|<7r. (C-4) 

This relation is important for understanding the analytic properties of the 
classical exponential integral function in that it isolates the multi-valued part 
represented by the logarithmic function from the regular part represented by 
the entire function Ein (z). In 1R + the function Ein (x) turns out to be a 
Bernstein function, which means that is positive, increasing, with the first 
derivative completely monotonic. 

The asymptotic behaviour as z — > 00 of the exponential integrals can be 
obtained from the integral representation (C.I) noticing that 



00 e 



-t 

-z 



EAz):= / dt = e~ z / du. (C.5) 

y 1 Jz t Jo u + z K J 

In fact, by repeated partial integrations in the R.H.S. , we get for |argz| < 

7T — 5 , 

e~ z 00 71 ! 

Z n=0 Z 

We now report a number of relevant Laplace transform pairs in which loga- 
rithmic and exponential integral functions are involved. 

Taking t > 0, the basic Laplace transforms pairs are 

£{log t; s} = - 7 + 1 ° gS , Ms > , (C.7) 

s 

C{£i(t);s} = log(g + 1) Sb>-1. (C.8) 
s 

The proof of (C.7) and (C.8) is found, for example, in the treatise by Ghizzetti 
and Ossicini [19], see Eqs. [4.6.15-16], pp. 104-105. We then easily derive 
for 9fo > 

£{ 7 + logt;s} = -^, (C.9) 
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rs ,i + , c u\ \ log(s + l) lo S s Iog(l/s + l) . . 

£{7 + logt + £i(t);s} = = , (CIO) 

s s s 

C{-< + log t + e'£ l(t ); s} = _ ]S5£ = . (C .ll) 

s — 1 s s(s — 1) 

We outline the different asymptotic behaviour of the three functions £i(t), 
Ein (t) and 7+logt+e* £i(t) , for small argument (t — ► + ) and large argument 
(t -> +oo). By using Eqs (C.2), (C.4) and (C.6), we have 



£ x {t)) ~ (C.12) 




Ein(i)~<^ (C.13) 



tlog(l/t), t-+0+, 
7 + logt + e*fi(t) ~ -( (C.14) 

logt , i — ► +oo . 

We note that all the above asymptotic representations can be obtained from 
the Laplace transforms of the corresponding functions by invoking the Taube- 
rian theory for regularly varying functions (power functions multiplied by 
slowly varying function^) , a topic not so well known which is adequately 
treated in the treatise on Probability by Feller [17], see Chapter XIII. 5. 

We conclude this subsection pointing out the Laplace transform pair 

£{v(t,a);s} = — ± , 3te > . (C.15) 

where 

poo t a+T 

v(t, a) := / dr, a > -1 . (C.16) 

jo 1 (a + t + 1) 

For details on this transcendental function the reader is referred to the third 
volume of the Handbook of the Bateman Project [16], see in Chapter XVIII 
(devoted to the Miscellaneous functions) §18.3, pp. 217-224. 



6 Definition: We call a (measurable) positive function a(y), defined in a right neigh- 
bourhood of zero, slowly varying at zero if a(y) > and a(cy)/a(y) — ► 1 with y — ► for 
every c > 0. We call a (measurable) positive function b(y), defined in a neighbourhood of 
infinity, slowly varying at infinity if b(cy)/b(y) — > 1 with y — > oo for every c > 0. Examples: 
(logy) 7 with 7 e H and exp (logy/log logy). 
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